clear all 
set more off

global root "ENTER ROOT FOLDER HERE"
global output "$root/Tables/"
cd "$root"

use "Data/Analysis/census.dta"

*** Merge in treatment assignment and district-level datasets
ren b_district dcode // district of birth code used for merging into other files
merge m:1 dcode using "Data/District-level/treatment_map.dta", gen(merge_treatment_map) // Merges in treatment assignment + admin units from 1967 census (same as existed in 1966)
merge m:1 district_1967 using "Data/District-level/taxation.dta", gen(merge_tax) // Merges in Jensen & Mkama (1968) and Lee (1965) district-level data
merge m:1 dcode using "Data/District-level/schools.dta", gen(merge_schools) // Merges in admin data on schools existing prior to reform

sum merge_* // check correct merge, should be 3 for all
drop merge_*

gen bw=yob-1966 // cohort bandwidth variable used to restrict cohorts used in main analysis

encode district_1967, gen(dcode_1967)
encode region_1967, gen(rcode_1967)
encode district_1967_g, gen(dcode_1967_g) // FEs used in main analysis defined at (pre-reform parent district * treatment) level

gen group = g_parent // baseline control districts specification: all districts sharing same pre-reform parent district
gen treat = group==1 & yob>=1966
replace treat = . if group==.
gen post = yob>=year_reform

encode dcode, gen(dcode_2012) // modern districts used in supplementary analyses (Table A4)
ren b_region rcode_2012 // modern regions used in supplementary analyses (Table A4)
egen dcode_1967_y = concat(yob dcode_1967_g) // district-cohort clustering used in supplementary analyses (Table A4)

* alternative treatment definitions (used in Table A4)
gen treat_urban = g_urban==1 & yob>=1966 // alternative control districts specification: towns/former towns as per 1967 census
replace treat_urban = . if g_urban==.
gen treat_all = g_parent==1 & yob>=1966 // alternative control districts specification: no restrictions on control districts

gen treat_full = 0 // full set of reform waves after initial reform (used in Table A11)
replace treat_full = 1 if yob>=year_reform // district-specific reform years (also used in Figure A10)

* generate lead terms for first stage *
forvalues i = 1(1)5 {
	gen lead_`i' = yob >= 1966-`i' & group==1
	replace lead_`i' = . if group==.	
	label variable lead_`i' "\hspace{0.1cm} \textit{Reform}\$_{t+`i'}\$"
}

label variable treat "\hspace{0.1cm} \textit{Reform}"
label variable group "\hspace{0.1cm} Born in treated district"
label variable post "\hspace{0.1cm} Born after reform"
label variable treat "\hspace{0.1cm} \textit{Reform}"
label variable treat_urban "\hspace{0.1cm} \textit{Reform}"
label variable treat_all "\hspace{0.1cm} \textit{Reform}"

label variable treat_full "\hspace{0.1cm} \textit{Reform}"
label define lab_t 0 "\hspace{0.1cm} \textit{Control}"
label define lab_t 1 "\hspace{0.1cm} \textit{Reform}", add
label values treat_full lab_t

label define lab_wave 1 "'66"
label define lab_wave 2 "'80s", add
label define lab_wave 3 "'09", add
label values wave lab_wave

* generate standardized district-level controls within the analysis sample (for use in Table 4)
foreach var of varlist jm_* lee_* sch_* {
	egen z_`var' = std(`var') if abs(bw) <= 10 & treat!=.
}

label variable z_jm_rate_payers_cap "Share paying tax"
label variable z_jm_tax_cap "Tax per capita"
label variable z_lee_diff_tax "\$(\tau^{\text{Max}} - \tau^{\text{Min}})\$"
label variable z_sch_primary_n_66 "Primary schools"
label variable z_sch_secondary_any_66 "Secondary school"
label variable z_jm_gdp_cap "GDP per capita"
label variable z_jm_pop_dens "Population density"

cd "$root"
savesome if abs(bw)<=10 & group!=. using "Data/Analysis/census_bw10.dta", replace // save restricted analysis file (for use in make_figures.R)

global fe = "yob dcode_1967_g" // FEs defined by cohort of birth and district*treatment status of locality of birth (as defined pre-reform)
global clust = "dcode_1967_g" // Clustering defined at district*treatment status of locality of birth (as defined pre-reform)
global ctrl = "c.pl_male c.pl_male#(i.dcode_1967_g)" // Flexibly control for locality of birth * gender differences


********************************************************************************
*** Descriptive statistics (Table A1) ******************************************
********************************************************************************

cd "$output"

local desc_vars "age pl_male pl_tanzanian pl_alive_father pl_alive_mother"
local fs_vars "b_cert group post treat"
local out_vars "ed_primary ed_secondary ed_uni s_nih s_pens_priv s_pens_pub"

estpost summarize `desc_vars' `fs_vars' `out_vars' if abs(bw) <= 10 & group!=., listwise
est store A
estpost summarize `desc_vars' `fs_vars' `out_vars' if abs(bw) <= 10 & group==1, listwise
est store B
estpost summarize `desc_vars' `fs_vars' `out_vars' if abs(bw) <= 10 & group==0, listwise
est store C

local over "& \multicolumn{2}{c}{Both} & \multicolumn{2}{c}{Treated} & \multicolumn{2}{c}{Control} \\"
local lines "\cmidrule(lr{0.5em}){2-3} \cmidrule(lr{0.5em}){4-5} \cmidrule(lr{0.5em}){6-7}"	
local names "& Mean & SD & Mean & SD & Mean & SD \\ \midrule"	

esttab A B C using "Table A1.tex", replace posthead("`over'" "`lines'" "`names'")  ///
	nonum cells("mean(fmt(2) label({})) sd(fmt(2) label({}))") ///
	label nomtitles b(3) se f nogaps booktabs sfmt(3 0) compress ///
	refcat(age "\textbf{A. Sample characteristics}" b_cert "\textbf{B. First stage variables}" ed_primary "\textbf{C. Outcome variables}", nolabel)


****************************************************************************************************************************************************************
*** First stage (Tables 1, A4, A6, A7, A11) ********************************************************************************************************************
****************************************************************************************************************************************************************	
	
********************************************************************************
*** First stage (Table 1) ******************************************************
********************************************************************************

cd "$output"

label variable treat "\hspace{0.1cm} \textit{Reform} ($\beta^{FS}$)"

local cols "fs_1_1 fs_1_2 fs_1_3 fs_2_1 fs_2_2 fs_2_3"

local t1 ""
local t2 "i.rcode_1967#c.yob" // Region of birth linear trends 
local t3 "i.dcode_1967_g#c.yob" // District of birth (as defined pre-reform) linear trends

local trends1 "None"
local trends2 "Region"
local trends3 "District"

estimates clear

foreach spec in 1 2 3 {
	
	* baseline specifications (varying use of unit-linear time trends)
	quietly reghdfe b_cert treat `t`spec'' $ctrl if abs(bw) <= 10, absorb( $fe ) cluster( $clust )
		estadd local Trends "`trends`spec''"
		sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		test treat
		estadd scalar fs = `r(F)'
		est sto fs_1_`spec'
		
	* specification with lead terms (varying use of unit-linear time trends)
	local leads "lead_1 lead_2 lead_3 lead_4 lead_5"
	quietly reghdfe b_cert treat  `leads' `t`spec'' $ctrl if abs(bw) <= 10, absorb( $fe ) cluster( $clust )
		estadd local Trends "`trends`spec''"
		sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		test treat
		estadd scalar fs = `r(F)'
		est sto fs_2_`spec'

}

* Main first stage table (Table 1)
esttab `cols' using "Table 1.tex", replace nomtitles b(2) se label f nogaps booktabs scalars("Trends Time trends" "fs F-statistic" "DV_Mean Outcome mean" "N Observations") sfmt(2 1 2 0) noobs compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels keep(t* l*) substitute(\_ _) order(treat lead_*)

* Supplementary version including control variable coefficients
esttab `cols' using "With controls/Table 1.tex", replace nomtitles b(2) se label f nogaps booktabs scalars("Trends Time trends" "fs F-statistic" "DV_Mean Outcome mean" "N Observations") sfmt(2 1 2 0) noobs compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels keep(t* l* pl_*) substitute(\_ _) order(treat lead_* pl_male)


********************************************************************************
*** First stage robustness (Tables A4, A6, A7, A11) ****************************
********************************************************************************

estimates clear

label variable treat "\hspace{0.1cm} \textit{Reform}"

*** Table A4: First stage (Robustness)
********************************************************************************

local cols "fs_1_1 fs_1_2 fs_1_3 fs_2_1 fs_2_2 fs_2_3"

*** Panel A (Varying included cohorts) *****
********************************************

local t1 ""
local t2 "i.rcode_1967#c.yob"
local t3 "i.dcode_1967_g#c.yob"

local trends1 "None"
local trends2 "Region"
local trends3 "District"

foreach spec in 1 2 3 {
	
	eststo fs_1_`spec': quietly reghdfe b_cert treat `t`spec'' $ctrl if abs(bw) <= 5, absorb( $fe ) cluster( $clust )
		estadd local Trends "`trends`spec''"
		quietly sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
	eststo fs_2_`spec': quietly reghdfe b_cert treat `t`spec'' $ctrl, absorb( $fe ) cluster( $clust )
		estadd local Trends "`trends`spec''"
		quietly sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'	

}

local titles "& \multicolumn{3}{c}{\textbf{+/- 5 cohorts}} & \multicolumn{3}{c}{\textbf{All cohorts}} \\"
local lines "\cmidrule(lr{0.5em}){2-4} \cmidrule(lr{0.5em}){5-7} "
local numbers "\textbf{A. Varying included cohorts} & (1) & (2) & (3) & (4) & (5) & (6) \\ \midrule"
esttab `cols' using "Table A4.tex", nonumber replace posthead("`titles'" "`lines'" "`numbers'") nomtitles b(2) se label f nogaps booktabs scalars("Trends Time trends" "DV_Mean Outcome mean" "N Observations") sfmt(3 2 0) noobs compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels keep(treat)
esttab `cols' using "With controls/Table A4.tex", nonumber replace posthead("`titles'" "`lines'" "`numbers'") nomtitles b(2) se label f nogaps booktabs scalars("Trends Time trends" "DV_Mean Outcome mean" "N Observations") sfmt(3 2 0) noobs compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels keep(treat pl_*)	
	
*** Panel B (Excluding birth years) ********
********************************************

foreach spec in 1 2 3 {
	
	eststo fs_1_`spec': reghdfe b_cert treat `t`spec'' $ctrl if abs(bw) <= 10 & yob != 1966, absorb( $fe ) cluster( $clust ) // excluding people born in reform year
		estadd local Trends "`trends`spec''"
		quietly sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
	eststo fs_2_`spec': reghdfe b_cert treat `t`spec'' $ctrl if abs(bw) <= 10 & age!=40&age!=45&age!=50&age!=55, absorb( $fe ) cluster( $clust ) // excluding people reporting "heaped" ages
		estadd local Trends "`trends`spec''"
		quietly sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'	

}

local titles "\midrule & \multicolumn{3}{c}{\textbf{-Reform year}} & \multicolumn{3}{c}{\textbf{-Heaped ages}} \\"
local numbers "\textbf{B. Excluding birth years} & (1) & (2) & (3) & (4) & (5) & (6) \\ \midrule"
esttab `cols' using "Table A4.tex", nonumber append posthead("`titles'" "`lines'" "`numbers'") nomtitles b(2) se label f nogaps booktabs scalars("Trends Time trends" "DV_Mean Outcome mean" "N Observations") sfmt(3 2 0) noobs compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels keep(treat)
esttab `cols' using "With controls/Table A4.tex", nonumber append posthead("`titles'" "`lines'" "`numbers'") nomtitles b(2) se label f nogaps booktabs scalars("Trends Time trends" "DV_Mean Outcome mean" "N Observations") sfmt(3 2 0) noobs compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels keep(treat pl_*)

*** Panel C (Different control variables) **
********************************************	

foreach spec in 1 2 3 {
	
	local distvars "c.z_jm_*"
	eststo fs_1_`spec': quietly reghdfe b_cert treat c.post##(`distvars') `t`spec'' $ctrl if abs(bw) <= 10, absorb( $fe ) cluster( $clust )
		estadd local Trends "`trends`spec''"
		quietly sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		
	local indvars "(c.pl_*) (c.pl_*)#i.dcode_1967_g"
	eststo fs_2_`spec': quietly reghdfe b_cert treat `indvars' `t`spec'' if abs(bw) <= 10, absorb( $fe ) cluster( $clust )
		estadd local Trends "`trends`spec''"
		quietly sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'	

}

local titles "\midrule & \multicolumn{3}{c}{\textbf{District-level}} & \multicolumn{3}{c}{\textbf{Individual-level}} \\"
local numbers "\textbf{C. Control variables} & (1) & (2) & (3) & (4) & (5) & (6) \\ \midrule"
esttab `cols' using "Table A4.tex", nonumber append posthead("`titles'" "`lines'" "`numbers'") nomtitles b(2) se label f nogaps booktabs scalars("Trends Time trends" "DV_Mean Outcome mean" "N Observations") sfmt(3 2 0) noobs compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels keep(treat)
esttab `cols' using "With controls/Table A4.tex", nonumber append posthead("`titles'" "`lines'" "`numbers'") nomtitles b(2) se label f nogaps booktabs scalars("Trends Time trends" "DV_Mean Outcome mean" "N Observations") sfmt(3 2 0) noobs compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels keep(treat pl_*)

*** Panel D (Different control districts) **
********************************************

foreach spec in 1 2 3 {
	
	cap drop treat_1
	gen treat_1 = treat_urban
	label variable treat_1 "\hspace{0.1cm} \textit{Reform}"

	eststo fs_1_`spec': quietly reghdfe b_cert treat_1 `t`spec'' $ctrl if abs(bw) <= 10, absorb( $fe ) cluster( $clust )
		estadd local Trends "`trends`spec''"
		quietly sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
	
	cap drop treat_1
	gen treat_1 = treat_all
	label variable treat_1 "\hspace{0.1cm} \textit{Reform}"	
		
	eststo fs_2_`spec': quietly reghdfe b_cert treat_1 `t`spec'' $ctrl if abs(bw) <= 10, absorb( $fe ) cluster( $clust )
		estadd local Trends "`trends`spec''"
		quietly sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'	

}
	
local titles "\midrule & \multicolumn{3}{c}{\textbf{Urban}} & \multicolumn{3}{c}{\textbf{Unrestricted}} \\"
local numbers "\textbf{D. Changing control districts} & (1) & (2) & (3) & (4) & (5) & (6) \\ \midrule"
esttab `cols' using "Table A4.tex", nonumber append posthead("`titles'" "`lines'" "`numbers'") nomtitles b(2) se label f nogaps booktabs scalars("Trends Time trends" "DV_Mean Outcome mean" "N Observations") sfmt(3 2 0) noobs compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels keep(treat*)
esttab `cols' using "With controls/Table A4.tex", nonumber append posthead("`titles'" "`lines'" "`numbers'") nomtitles b(2) se label f nogaps booktabs scalars("Trends Time trends" "DV_Mean Outcome mean" "N Observations") sfmt(3 2 0) noobs compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels keep(treat* pl_*)

*** Panel E (Different estimation using modern FEs or more granular clustering) **
********************************************

foreach spec in 1 2 3 {

	local t1 ""
	local t2 "i.rcode_2012#c.yob"
	local t3 "i.dcode_2012#c.yob"

	eststo fs_1_`spec': quietly reghdfe b_cert treat `t`spec'' $ctrl if abs(bw) <= 10, absorb( yob dcode_2012 ) cluster( dcode_2012 )
		estadd local Trends "`trends`spec''"
		quietly sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
	eststo fs_2_`spec': quietly reghdfe b_cert treat `t`spec'' $ctrl if abs(bw) <= 10, absorb( $fe ) cluster( dcode_1967_y )
		estadd local Trends "`trends`spec''"
		quietly sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'	

}
	
local titles "\midrule & \multicolumn{3}{c}{\textbf{2012 district FEs}} & \multicolumn{3}{c}{\textbf{District-cohort clustering}} \\"
local numbers "\textbf{E. Varying estimation} & (1) & (2) & (3) & (4) & (5) & (6) \\ \midrule"
esttab `cols' using "Table A4.tex", nonumber append posthead("`titles'" "`lines'" "`numbers'") nomtitles b(2) se label f nogaps booktabs scalars("Trends Time trends" "DV_Mean Outcome mean" "N Observations") sfmt(3 2 0) noobs compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels keep(treat*)
esttab `cols' using "With controls/Table A4.tex", nonumber append posthead("`titles'" "`lines'" "`numbers'") nomtitles b(2) se label f nogaps booktabs scalars("Trends Time trends" "DV_Mean Outcome mean" "N Observations") sfmt(3 2 0) noobs compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels keep(treat* pl_*)


*** Panel F (NPS Dataset - see make_tables_nps.do) ******
********************************************

	

*** Table A6: First stage (Alternative estimation)
********************************************************************************

local cols "fs_2 fs_5 fs_10 fs_15 fs_20 fs_30"

*** Panel A (Regression discontinuity-style specification within treated districts)
********************************************

foreach cohorts in 2 5 10 15 20 30 {
	eststo fs_`cohorts': quietly reghdfe b_cert c.post##c.bw $ctrl if abs(bw)<=`cohorts' & group==1, a( $clust ) cluster( $clust )
		estadd local BW "`cohorts'"	
		sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'	
}	

local numbers "\textbf{A. Regression discontinuity} & (1) & (2) & (3) & (4) & (5) & (6) \\ \midrule"
esttab `cols' using "Table A6.tex", nonumber replace posthead("`numbers'") nomtitles b(2) se label f nogaps booktabs scalars("BW Bandwidth" "DV_Mean Outcome mean" "N Observations") sfmt(0 2 0) noobs compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels keep(post)
esttab `cols' using "With controls/Table A6.tex", nonumber replace posthead("`numbers'") nomtitles b(2) se label f nogaps booktabs scalars("BW Bandwidth" "DV_Mean Outcome mean" "N Observations") sfmt(0 2 0) noobs compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels keep(post pl_*)


*** Panel B (Household FEs) ****************
********************************************

foreach cohorts in 2 5 10 15 20 30 {
	eststo fs_`cohorts': quietly reghdfe b_cert treat $ctrl if abs(bw)<=`cohorts', a(i_hhid $fe ) cluster( $clust )
		estadd local BW "`cohorts'"	
		sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'	
}
	
local numbers "\midrule \textbf{B. Household fixed effects} & (1) & (2) & (3) & (4) & (5) & (6) \\ \midrule"
esttab `cols' using "Table A6.tex", nonumber append posthead("`numbers'") nomtitles b(2) se label f nogaps booktabs scalars("BW Bandwidth" "DV_Mean Outcome mean" "N Observations") sfmt(0 2 0) noobs compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels keep(treat)
esttab `cols' using "With controls/Table A6.tex", nonumber append posthead("`numbers'") nomtitles b(2) se label f nogaps booktabs scalars("BW Bandwidth" "DV_Mean Outcome mean" "N Observations") sfmt(0 2 0) noobs compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels keep(treat pl_*)

	
*** Table A7: Placebo outcomes
********************************************************************************

foreach var in pl_male pl_tanzanian pl_alive_father pl_alive_mother {
	eststo `var': capture quietly reghdfe `var' treat if abs(bw) <= 10, absorb( $fe ) cluster( $clust )
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
}

local vars "& Male & Tanzanian & \specialcell{Father \\ alive} & \specialcell{Mother \\ alive} \\"
local numbers " & (1) & (2) & (3) & (4) \\ \midrule"

esttab pl_male pl_tanzanian pl_alive_father pl_alive_mother using "Table A7.tex", ///
	replace nonumber posthead("`vars'" "`numbers'") nomtitles b(2) se label f nogaps booktabs scalars("DV_Mean Outcome mean" "N Observations") sfmt(2 0) noobs compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels noconstant keep(treat)
	
	
*** Table A11: Comparison of registration reforms (use modern districts to avoid within-FE variation in reform date)
********************************************************************************

foreach var in b_cert {	
	
	eststo `var'_1: capture quietly reghdfe `var' 1.treat_full##i.wave $ctrl, absorb(yob dcode_2012) cluster(dcode_2012) // pooled comparison of reforms
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
	eststo `var'_2: capture quietly reghdfe `var' 1.treat_full $ctrl if abs(bw)<=10, absorb(yob dcode_2012) cluster(dcode_2012) // baseline first stage (all control districts)
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
	eststo `var'_3: capture quietly reghdfe `var' 1.treat_full $ctrl if yob>=1971 & yob<1998, absorb(yob dcode_2012) cluster(dcode_2012) // 1980s reform wave
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
	eststo `var'_4: capture quietly reghdfe `var' 1.treat_full $ctrl if yob>=1999, absorb(yob dcode_2012) cluster(dcode_2012) // 2000s reform wave
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
}

local yrs " & Pooled & '66 & '80s & '09 \\ \midrule"
local numbers " & (1) & (2) & (3) & (4) \\ \midrule"
esttab b_cert_1 b_cert_2 b_cert_3 b_cert_4 using "Table A11.tex", ///
	replace nonumber posthead("`yrs'" "`numbers'") nomtitles b(2) se label f nogaps booktabs scalars("DV_Mean Outcome mean" "N Observations") sfmt(2 0) noobs compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels keep(*at_ful*) substitute(\_ _)		
esttab b_cert_1 b_cert_2 b_cert_3 b_cert_4 using "With controls/Table A11.tex", ///
	replace nonumber posthead("`yrs'" "`numbers'") nomtitles b(2) se label f nogaps booktabs scalars("DV_Mean Outcome mean" "N Observations") sfmt(2 0) noobs compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels keep(*at_ful* pl_*) substitute(\_ _)		
	
	

****************************************************************************************************************************************************************
*** IV results (Table 2, A8, A9, A10) **************************************************************************************************************************
****************************************************************************************************************************************************************

keep if abs(bw) <= 10 // Keep only baseline sample of cohorts +/- 10 years. 
* Comment above line out if trying to estimate results with no restrictions on cohorts or larger set of cohorts (i.e. additional specification A2 below) but that estimation will take a long time

local ed_vars "ed_primary ed_secondary ed_uni" // used in Table 2
local lit_vars "lit_any lit_kiswa lit_eng" // used in Table A8
local ss_vars "s_nih s_pens_priv s_pens_pub"  // used in Table 2
local s_vars "s_any s_nih s_nssf s_pension_parastatal s_pension_pubservice s_pension_govemp s_pension_locauth s_other" // used in Table A9

foreach var in `ed_vars' `lit_vars' `ss_vars' `s_vars' {
	capture drop z_`var'
	egen z_`var' = std(`var') if abs(bw) <= 10 
}

foreach var in `ed_vars' `ss_vars' `lit_vars' `s_vars' { 
	
	eststo `var': quietly reghdfe `var' b_cert $ctrl if abs(bw)<=10 & group!=., absorb( $fe ) cluster( $clust )	// OLS estimation
	quietly ivreghdfe `var' (b_cert=treat) $ctrl i.yob if abs(bw)<=10 & group!=., a( $clust ) cluster( $clust ) // IV estimation
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		matrix define pmat = J(1,2,0)
		matrix pmat[1,2] = 2*ttail(e(df_r),abs(_b[pl_male]/_se[pl_male]))
		matrix colnames pmat = b_cert pl_male
		boottest b_cert, reps(5000) seed(94305) noci statistic(t)
		matrix pmat[1,1] = `r(p)'
		estadd matrix pboot = pmat
		est sto `var'_iv
	eststo z_`var': quietly reghdfe z_`var' z_asset_pca $ctrl if abs(bw)<=10 & group!=., absorb( $fe ) cluster( $clust ) // rho vector
		
}

cd "$output"

*** Table 2: Effects on access to the state
********************************************************************************

local titles "& \multicolumn{3}{c}{I. Education} & \multicolumn{3}{c}{II. Social security} \\"
local lines "\cmidrule(lr{0.5em}){2-4} \cmidrule(lr{0.5em}){5-7} "
local names " & Pri. & Sec. & Uni. & HI & Priv. & State \\"
local numbers " & (1) & (2) & (3) & (4) & (5) & (6) \\ \midrule"

esttab ed_primary ed_secondary ed_uni s_nih s_pens_priv s_pens_pub using "Table 2.tex", coeflabel(b_cert "\hspace{0.1cm} Registered ($\beta^{OLS})$") replace nonumber posthead("`titles'" "`lines'" "`names'" "`numbers'") nomtitles b(2) se label f nogaps booktabs sfmt(0) noobs compress nobaselevels keep(b_cert) star(* 0.10 ** 0.05 *** 0.01)
esttab ed_primary_iv ed_secondary_iv ed_uni_iv s_nih_iv s_pens_priv_iv s_pens_pub_iv using "Table 2.tex", append cells(b(star pvalue(pboot) fmt(%9.2f)) se(par fmt(%9.2f)) pboot(par([ ]) fmt(%9.2f))) coeflabel(b_cert  "\hspace{0.1cm} $\widehat{\text{Registered}}$ ($\beta^{IV})$") nonumber nomtitles label f nogaps booktabs noobs stats(DV_Mean widstat, labels("\midrule DV Mean" "FS F-statistic") layout("@") fmt(2 1)) compress nobaselevels keep(b_cert) nolines prehead(\midrule) star(* 0.10 ** 0.05 *** 0.01) mlabels(,none) collabels(,none)
esttab z_ed_primary z_ed_secondary z_ed_uni z_s_nih z_s_pens_priv z_s_pens_pub using "Table 2.tex", coeflabel(z_asset_pca "\$\rho\$(Wealth, DV)") append nonumber nomtitles b(2) not label f nogaps booktabs sfmt(0) noobs scalars("N Observations") compress nostar nobaselevels keep(z_asset_pca)

* Supplementary version including control variable coefficients
esttab ed_primary ed_secondary ed_uni s_nih s_pens_priv s_pens_pub using "With controls/Table 2.tex", coeflabel(b_cert "\hspace{0.1cm} Registered ($\beta^{OLS})$") replace nonumber posthead("`titles'" "`lines'" "`names'" "`numbers'") nomtitles b(2) se label f nogaps booktabs sfmt(0) noobs compress nobaselevels keep(b_cert pl_*) order(pl_* b_cert ) star(* 0.10 ** 0.05 *** 0.01)
esttab ed_primary_iv ed_secondary_iv ed_uni_iv s_nih_iv s_pens_priv_iv s_pens_pub_iv using "With controls/Table 2.tex", cells(b(star pvalue(pboot) fmt(%9.2f)) se(par fmt(%9.2f)) pboot(par([ ]) fmt(%9.2f))) coeflabel(b_cert "\hspace{0.1cm} $\widehat{\text{Registered}}$ ($\beta^{IV})$") append nonumber nomtitles label f nogaps booktabs noobs stats(DV_Mean widstat, labels("\midrule DV Mean" "FS F-statistic") layout("@") fmt(2 1)) compress nobaselevels keep(b_cert pl_*) order(pl_*  b_cert ) nolines prehead(\midrule) star(* 0.10 ** 0.05 *** 0.01) mlabels(,none) collabels(,none)
esttab z_ed_primary z_ed_secondary z_ed_uni z_s_nih z_s_pens_priv z_s_pens_pub using "With controls/Table 2.tex", coeflabel(z_asset_pca "\$\rho\$(Wealth, DV)") append nonumber nomtitles b(2) not label f nogaps booktabs sfmt(0) noobs scalars("N Observations") compress nostar nobaselevels keep(z_asset_pca)
	
	
*** Table A8: Effects on literacy
********************************************************************************

local names " & Any & Kisw. & Eng. \\"
local numbers " & (1) & (2) & (3) \\ \midrule"

esttab lit_any lit_kiswa lit_eng using "Table A8.tex", coeflabel(b_cert "\hspace{0.1cm} Registered ($\beta^{OLS})$") replace nonumber posthead("`names'" "`numbers'") nomtitles b(2) se label f nogaps booktabs sfmt(0) noobs compress nobaselevels keep(b_cert)	star(* 0.10 ** 0.05 *** 0.01)
esttab lit_any_iv lit_kiswa_iv lit_eng_iv using "Table A8.tex", append cells(b(star pvalue(pboot) fmt(%9.2f)) se(par fmt(%9.2f)) pboot(par([ ]) fmt(%9.2f))) coeflabel(b_cert  "\hspace{0.1cm} $\widehat{\text{Registered}}$ ($\beta^{IV})$") nonumber nomtitles label f nogaps booktabs sfmt(2 1 0) noobs stats(DV_Mean widstat, labels("\midrule DV Mean" "FS F-statistic") layout("@") fmt(2 2 1)) compress nobaselevels keep(b_cert) nolines prehead(\midrule) star(* 0.10 ** 0.05 *** 0.01) mlabels(,none) collabels(,none)
esttab z_lit_any z_lit_kiswa z_lit_eng using "Table A8.tex", coeflabel(z_asset_pca "\$\rho\$(Wealth, DV)") append nonumber nomtitles b(2) not label f nogaps booktabs sfmt(0) noobs scalars("N Observations") compress nostar nobaselevels keep(z_asset_pca)

* Supplementary version including control variable coefficients
esttab lit_any lit_kiswa lit_eng using "With controls/Table A8.tex", coeflabel(b_cert "\hspace{0.1cm} Registered ($\beta^{OLS})$") replace nonumber posthead("`names'" "`numbers'") nomtitles b(2) se label f nogaps booktabs sfmt(0) noobs compress nobaselevels keep(b_cert pl_*) order(pl_* b_cert) star(* 0.10 ** 0.05 *** 0.01)
esttab lit_any_iv lit_kiswa_iv lit_eng_iv using "With controls/Table A8.tex", cells(b(star pvalue(pboot) fmt(%9.2f)) se(par fmt(%9.2f)) pboot(par([ ]) fmt(%9.2f))) coeflabel(b_cert "\hspace{0.1cm} $\widehat{\text{Registered}}$ ($\beta^{IV})$") append nonumber nomtitles label f nogaps booktabs noobs stats(DV_Mean widstat, labels("\midrule DV Mean" "FS F-statistic") layout("@") fmt(2 1)) compress nobaselevels keep(b_cert pl_*) order(pl_*  b_cert ) nolines prehead(\midrule) star(* 0.10 ** 0.05 *** 0.01) mlabels(,none) collabels(,none)
esttab z_lit_any z_lit_kiswa z_lit_eng using "With controls/Table A8.tex", coeflabel(z_asset_pca "\$\rho\$(Wealth, DV)") append nonumber nomtitles b(2) not label f nogaps booktabs sfmt(0) noobs scalars("N Observations") compress nostar nobaselevels keep(z_asset_pca)

*** Table A9: Effects on access to social security
********************************************************************************
	
local titles "& & & \multicolumn{2}{c}{Private} & \multicolumn{3}{c}{State} & \\"
local lines "\cmidrule(lr{0.5em}){4-5} \cmidrule(lr{0.5em}){6-8} "
local names " & Any & HI & NSSF & PPF & PSPF & GEPF & LAPF & Other \\"
local numbers " & (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8) \\ \midrule"

esttab s_any s_nih s_nssf s_pension_parastatal s_pension_pubservice s_pension_govemp s_pension_locauth s_other using "Table A9.tex", coeflabel(b_cert "\hspace{0.1cm} Registered ($\beta^{OLS})$") replace nonumber posthead("`titles'" "`lines'" "`names'" "`numbers'") nomtitles b(2) se label f nogaps booktabs sfmt(0) noobs compress nobaselevels keep(b_cert)	star(* 0.10 ** 0.05 *** 0.01)
esttab s_any_iv s_nih_iv s_nssf_iv s_pension_parastatal_iv s_pension_pubservice_iv s_pension_govemp_iv s_pension_locauth_iv s_other_iv using "Table A9.tex", append cells(b(star pvalue(pboot) fmt(%9.2f)) se(par fmt(%9.2f)) pboot(par([ ]) fmt(%9.2f))) coeflabel(b_cert  "\hspace{0.1cm} $\widehat{\text{Registered}}$ ($\beta^{IV})$") nonumber nomtitles label f nogaps booktabs noobs stats(DV_Mean widstat, labels("\midrule DV Mean" "FS F-statistic") layout("@") fmt(2 1)) compress nobaselevels keep(b_cert) nolines prehead(\midrule) star(* 0.10 ** 0.05 *** 0.01) mlabels(,none) collabels(,none)
esttab z_s_any z_s_nih z_s_nssf z_s_pension_parastatal z_s_pension_pubservice z_s_pension_govemp z_s_pension_locauth z_s_other using "Table A9.tex", coeflabel(z_asset_pca "\$\rho\$(Wealth, DV)") ///
	append nonumber nomtitles b(2) not label f nogaps booktabs sfmt(0) noobs scalars("N Observations") compress nostar nobaselevels keep(z_asset_pca)
	
* Supplementary version including control variable coefficients
esttab s_any s_nih s_nssf s_pension_parastatal s_pension_pubservice s_pension_govemp s_pension_locauth s_other using "With controls/Table A9.tex", coeflabel(b_cert "\hspace{0.1cm} Registered ($\beta^{OLS})$") replace nonumber posthead("`titles'" "`lines'" "`names'" "`numbers'") nomtitles b(2) se label f nogaps booktabs sfmt(0) noobs compress nobaselevels keep(b_cert pl_*) order(pl_* b_cert) star(* 0.10 ** 0.05 *** 0.01)
esttab s_any_iv s_nih_iv s_nssf_iv s_pension_parastatal_iv s_pension_pubservice_iv s_pension_govemp_iv s_pension_locauth_iv s_other_iv using "With controls/Table A9.tex", cells(b(star pvalue(pboot) fmt(%9.2f)) se(par fmt(%9.2f)) pboot(par([ ]) fmt(%9.2f))) coeflabel(b_cert "\hspace{0.1cm} $\widehat{\text{Registered}}$ ($\beta^{IV})$") append nonumber nomtitles label f nogaps booktabs noobs stats(DV_Mean widstat, labels("\midrule DV Mean" "FS F-statistic") layout("@") fmt(2 1)) compress nobaselevels keep(b_cert pl_*) order(pl_*  b_cert ) nolines prehead(\midrule) star(* 0.10 ** 0.05 *** 0.01) mlabels(,none) collabels(,none)
esttab z_s_any z_s_nih z_s_nssf z_s_pension_parastatal z_s_pension_pubservice z_s_pension_govemp z_s_pension_locauth z_s_other using "With controls/Table A9.tex", coeflabel(z_asset_pca "\$\rho\$(Wealth, DV)") ///
	append nonumber nomtitles b(2) not label f nogaps booktabs sfmt(0) noobs scalars("N Observations") compress nostar nobaselevels keep(z_asset_pca)
	
	
****************************************************************************************************************************************************************
*** Heterogeneity in compliance (Table 4) ********************************************************************************************************************
****************************************************************************************************************************************************************	

estimates clear

cd "$output"

foreach var in b_cert {

	eststo `var': quietly reghdfe `var' treat $ctrl if abs(bw) <= 10 & group!=., a( $fe ) cluster( $clust )
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
	
	* Panel I: Taxation
	eststo `var'_1_1: quietly reghdfe `var' c.treat##(c.z_jm_rate_payers_cap) $ctrl c.treat#($het) if abs(bw) <= 10 & group!=., a( $fe ) cluster( $clust )
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
	eststo `var'_1_2: quietly reghdfe `var' c.treat##(c.z_jm_tax_cap) $ctrl c.treat#($het) if abs(bw) <= 10 & group!=., a( $fe ) cluster( $clust )
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
	eststo `var'_1_3: quietly reghdfe `var' c.treat##(c.z_lee_diff_tax) $ctrl c.treat#($het) if abs(bw) <= 10 & group!=., a( $fe ) cluster( $clust )
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
	eststo `var'_1_4: quietly reghdfe `var' c.treat##(c.z_jm_rate_payers_cap c.z_jm_tax_cap c.z_lee_diff_tax) $ctrl c.treat#($het) if abs(bw) <= 10 & group!=., a( $fe ) cluster( $clust )
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
	
	* Panel II: Local public goods
	eststo `var'_2_1: quietly reghdfe `var' c.treat##(c.z_sch_primary_n_66) $ctrl c.treat#($het) if abs(bw) <= 10 & group!=., a( $fe ) cluster( $clust )
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
	eststo `var'_2_2: quietly reghdfe `var' c.treat##(c.z_sch_secondary_any_66) $ctrl c.treat#($het) if abs(bw) <= 10 & group!=., a( $fe ) cluster( $clust )
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
	eststo `var'_2_3: quietly reghdfe `var' c.treat##(c.z_sch_primary_n_66 c.z_sch_secondary_any_66) $ctrl c.treat#($het) if abs(bw) <= 10 & group!=., a( $fe ) cluster( $clust )
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
	
}

local over "& & \multicolumn{4}{c}{I. Taxation} & \multicolumn{3}{c}{II. Local public goods} \\"
local lines "\cmidrule(lr{0.5em}){3-6} \cmidrule(lr{0.5em}){7-9}"
local numbers "& (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8) \\ \midrule"
esttab b_cert b_cert_1_1 b_cert_1_2 b_cert_1_3 b_cert_1_4 b_cert_2_1 b_cert_2_2 b_cert_2_3 using "Table 4.tex", ///
	replace posthead("`over'" "`lines'" "`numbers'") nomtitles b(3) se label f nogaps booktabs scalars("DV_Mean Outcome mean" "N Observations") /// 
	sfmt(2 0) noobs nonumber compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels  substitute(\_ _) /// 
	keep(treat *treat#c.z_*) drop(*treat#c.z_jm_gdp_* *treat#c.z_jm_pop_*)

esttab b_cert b_cert_1_1 b_cert_1_2 b_cert_1_3 b_cert_1_4 b_cert_2_1 b_cert_2_2 b_cert_2_3 using "With controls/Table 4.tex", ///
	replace posthead("`over'" "`lines'" "`numbers'") nomtitles b(3) se label f nogaps booktabs scalars("DV_Mean Outcome mean" "N Observations") /// 
	sfmt(2 0) noobs nonumber compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels  substitute(\_ _) /// 
	keep(treat *treat#c.z_* pl_male) order(treat *treat#c.z_jm_rate* *treat#c.z_jm_tax* *treat#z_lee* treat#z_sch* *treat#c.z_jm_gdp_* *treat#c.z_jm_pop_* pl_male)
	
	
****************************************************************************************************************************************************************
*** Complier characteristics (Table 5) *************************************************************************************************************************
****************************************************************************************************************************************************************	

* See make_figures.R *
	
	
****************************************************************************************************************************************************************
*** Supplementary results (reported in additional_results.pdf) *************************************************************************************************
****************************************************************************************************************************************************************

estimates clear

*** Table 1 using every estimating permutation from Table A4 
********************************************************************************

forvalues i = 1(1)5 { // generating lead terms for both treat_urban and treat_all robustness tests (specification D)
	gen lead1_`i' = yob >= 1966-`i' & group==1
	replace lead1_`i' = . if g_urban==.	
	label variable lead1_`i' "\hspace{0.1cm} \textit{Reform}\$_{t+`i'}\$"
	
	gen lead2_`i' = yob >= 1966-`i' & group==1
	label variable lead2_`i' "\hspace{0.1cm} \textit{Reform}\$_{t+`i'}\$"
}

local t1 ""
local t2 "i.rcode_1967#c.yob"
local t3 "i.dcode_1967_g#c.yob"
local t4 "lead_1 lead_2 lead_3 lead_4 lead_5"
local t5 "lead_1 lead_2 lead_3 lead_4 lead_5 i.rcode_1967#c.yob"
local t6 "lead_1 lead_2 lead_3 lead_4 lead_5 i.dcode_1967_g#c.yob"

local trends1 "None"
local trends2 "Region"
local trends3 "District"
local trends4 "None"
local trends5 "Region"
local trends6 "District"

foreach spec in 1 2 3 4 5 6 {
	
	eststo A_1_`spec': quietly reghdfe b_cert treat `t`spec'' $ctrl if abs(bw) <= 5 & group!=., absorb( $fe ) cluster( $clust )
		estadd local Trends "`trends`spec''"
		quietly sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		test treat
		estadd scalar fs = `r(F)'
	eststo A_2_`spec': quietly reghdfe b_cert treat `t`spec'' $ctrl if group!=., absorb( $fe ) cluster( $clust )
		estadd local Trends "`trends`spec''"
		quietly sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'	
		test treat
		estadd scalar fs = `r(F)'
		
	eststo B_1_`spec': quietly reghdfe b_cert treat `t`spec'' $ctrl if bw!=0 & abs(bw)<=10 & group!=., absorb( $fe ) cluster( $clust )
		estadd local Trends "`trends`spec''"
		quietly sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		test treat
		estadd scalar fs = `r(F)'
	eststo B_2_`spec': quietly reghdfe b_cert treat `t`spec'' $ctrl if age!=40&age!=45&age!=50&age!=55 & abs(bw)<=10 & group!=., absorb( $fe ) cluster( $clust )
		estadd local Trends "`trends`spec''"
		quietly sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'	
		test treat
		estadd scalar fs = `r(F)'
	
	local distvars "c.z_jm_*"
	eststo C_1_`spec': quietly reghdfe b_cert treat `t`spec'' c.post#(`distvars') $ctrl if abs(bw)<=10 & group!=., absorb( $fe ) cluster( $clust )
		estadd local Trends "`trends`spec''"
		quietly sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		test treat
		estadd scalar fs = `r(F)'
	local indvars "c.pl_* i.dcode_2012#(c.pl_*)"
	eststo C_2_`spec': quietly reghdfe b_cert treat `t`spec'' `indvars' $ctrl if abs(bw)<=10 & group!=., absorb( $fe ) cluster( $clust )
		estadd local Trends "`trends`spec''"
		quietly sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		test treat
		estadd scalar fs = `r(F)'
	
	local t4 "lead1_1 lead1_2 lead1_3 lead1_4 lead1_5"
	local t5 "lead1_1 lead1_2 lead1_3 lead1_4 lead1_5 i.rcode_1967#c.yob"
	local t6 "lead1_1 lead1_2 lead1_3 lead1_4 lead1_5 i.dcode_1967_g#c.yob"
	eststo D_1_`spec': quietly reghdfe b_cert treat_urban `t`spec'' $ctrl if abs(bw)<=10 & g_urban!=., absorb( $fe ) cluster( $clust )
		estadd local Trends "`trends`spec''"
		quietly sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		test treat_urban
		estadd scalar fs = `r(F)'
	local t4 "lead2_1 lead2_2 lead2_3 lead2_4 lead2_5"
	local t5 "lead2_1 lead2_2 lead2_3 lead2_4 lead2_5 i.rcode_1967#c.yob"
	local t6 "lead2_1 lead2_2 lead2_3 lead2_4 lead2_5 i.dcode_1967_g#c.yob"
	eststo D_2_`spec': quietly reghdfe b_cert treat_all `t`spec'' $ctrl if abs(bw)<=10, absorb( $fe ) cluster( $clust )
		estadd local Trends "`trends`spec''"
		quietly sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		test treat_all
		estadd scalar fs = `r(F)'
			
	local t2 "i.rcode_2012#c.yob"
	local t3 "i.dcode_2012#c.yob"
	local t4 "lead_1 lead_2 lead_3 lead_4 lead_5"
	local t5 "lead_1 lead_2 lead_3 lead_4 lead_5 i.rcode_1967#c.yob"
	local t6 "lead_1 lead_2 lead_3 lead_4 lead_5 i.dcode_1967_g#c.yob"
	eststo E_1_`spec': quietly reghdfe b_cert treat `t`spec'' $ctrl if abs(bw)<=10 & group!=., absorb( yob dcode_2012 ) cluster( dcode_2012 )
		estadd local Trends "`trends`spec''"
		quietly sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		test treat
		estadd scalar fs = `r(F)'
	local t2 "i.rcode_1967#c.yob"
	local t3 "i.dcode_1967_g#c.yob"		
	eststo E_2_`spec': quietly reghdfe b_cert treat `t`spec'' $ctrl if abs(bw)<=10 & group!=., absorb( $fe ) cluster( dcode_1967_y )
		estadd local Trends "`trends`spec''"
		quietly sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		test treat
		estadd scalar fs = `r(F)'

}

foreach v1 in A B C D E {
	foreach v2 in 1 2 {
		esttab `v1'_`v2'_1 `v1'_`v2'_2 `v1'_`v2'_3 `v1'_`v2'_4 `v1'_`v2'_5 `v1'_`v2'_6 using "Additional results/Table 1_`v1'`v2'.tex", replace nomtitles b(2) se label f nogaps booktabs scalars("Trends Time trends" "fs F-statistic" "DV_Mean Outcome mean" "N Observations") sfmt(2 1 2 0) noobs compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels keep(t* l*) substitute(\_ _) order(t* le*)
	}
}

*** Table 2 using every estimating permutation from Table A4 
********************************************************************************

local ed_vars "ed_primary ed_secondary ed_uni"
local ss_vars "s_nih s_pens_priv s_pens_pub"

foreach var in `ed_vars' `ss_vars' {
	
	// A: Varying cohorts
	eststo A_1_`var': capture quietly reghdfe `var' b_cert $ctrl if abs(bw)<=5 & group!=., absorb( $fe ) cluster( $clust )	
	quietly ivreghdfe `var' (b_cert=treat) $ctrl i.yob if abs(bw)<=5 & group!=., a( $clust ) cluster( $clust )
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		matrix define bpmat = J(1,1,0)
		matrix colnames bpmat = "b_cert"
		boottest b_cert, reps(5000) seed(94305) noci statistic(t)
		matrix bpmat[1,1] = `r(p)'
		estadd matrix pboot = bpmat
		est sto A_1_`var'_iv
	eststo A_1_z_`var': quietly reghdfe z_`var' z_asset_pca $ctrl if abs(bw)<=5 & group!=., absorb( $fe ) cluster( $clust )
	
	eststo A_2_`var': capture quietly reghdfe `var' b_cert $ctrl if group!=., absorb( $fe ) cluster( $clust )	
	quietly ivreghdfe `var' (b_cert=treat) $ctrl i.yob if group!=., a( $clust ) cluster( $clust )
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		matrix define bpmat = J(1,1,0)
		matrix colnames bpmat = "b_cert"
		boottest b_cert, reps(5000) seed(94305) noci statistic(t)
		matrix bpmat[1,1] = `r(p)'
		estadd matrix pboot = bpmat
		est sto A_2_`var'_iv
	eststo A_2_z_`var': quietly reghdfe z_`var' z_asset_pca $ctrl if group!=., absorb( $fe ) cluster( $clust )
	
	// B: Excluding birth years
	eststo B_1_`var': capture quietly reghdfe `var' b_cert $ctrl if abs(bw)<=10 & bw!=0 & group!=., absorb( $fe ) cluster( $clust )	
	quietly ivreghdfe `var' (b_cert=treat) $ctrl i.yob if abs(bw)<=10 & bw!=0 & group!=., a( $clust ) cluster( $clust )
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		matrix define bpmat = J(1,1,0)
		matrix colnames bpmat = "b_cert"
		boottest b_cert, reps(5000) seed(94305) noci statistic(t)
		matrix bpmat[1,1] = `r(p)'
		estadd matrix pboot = bpmat
		est sto B_1_`var'_iv	
	eststo B_1_z_`var': quietly reghdfe z_`var' z_asset_pca $ctrl if abs(bw)<=10 & bw!=0 & group!=., absorb( $fe ) cluster( $clust )
	
	eststo B_2_`var': capture quietly reghdfe `var' b_cert $ctrl if abs(bw)<=10 & group!=. & age!=40&age!=45&age!=50&age!=55, absorb( $fe ) cluster( $clust )	
	quietly ivreghdfe `var' (b_cert=treat) $ctrl i.yob if abs(bw)<=10 & group!=. & age!=40&age!=45&age!=50&age!=55, a( $clust ) cluster( $clust )
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		matrix define bpmat = J(1,1,0)
		matrix colnames bpmat = "b_cert"
		boottest b_cert, reps(5000) seed(94305) noci statistic(t)
		matrix bpmat[1,1] = `r(p)'
		estadd matrix pboot = bpmat
		est sto B_2_`var'_iv		
	eststo B_2_z_`var': quietly reghdfe z_`var' z_asset_pca $ctrl if abs(bw)<=10 & group!=. & age!=40&age!=45&age!=50&age!=55, absorb( $fe ) cluster( $clust )
		
	// C: Control variables
	local distvars "c.z_jm_*"
	eststo C_1_`var': capture quietly reghdfe `var' b_cert c.post#(`distvars') $ctrl if abs(bw)<=10 & group!=., absorb( $fe ) cluster( $clust )	
	quietly ivreghdfe `var' (b_cert=treat) c.post#(`distvars') $ctrl i.yob if abs(bw)<=10 & group!=., a( $clust ) cluster( $clust )
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		matrix define bpmat = J(1,1,0)
		matrix colnames bpmat = "b_cert"
		boottest b_cert, reps(5000) seed(94305) noci statistic(t)
		matrix bpmat[1,1] = `r(p)'
		estadd matrix pboot = bpmat
		est sto C_1_`var'_iv
	eststo C_1_z_`var': quietly reghdfe z_`var' z_asset_pca c.post#(`distvars') $ctrl if abs(bw)<=10 & group!=., absorb( $fe ) cluster( $clust )
	
	local indvars "c.pl_* i.dcode_2012#(c.pl_*)"
	eststo C_2_`var': capture quietly reghdfe `var' b_cert `indvars' if abs(bw)<=10 & group!=., absorb( $fe ) cluster( $clust )	
	quietly ivreghdfe `var' (b_cert=treat) `indvars' i.yob if abs(bw)<=10 & group!=., a( $clust ) cluster( $clust )
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		matrix define bpmat = J(1,1,0)
		matrix colnames bpmat = "b_cert"
		boottest b_cert, reps(5000) seed(94305) noci statistic(t)
		matrix bpmat[1,1] = `r(p)'
		estadd matrix pboot = bpmat
		est sto C_2_`var'_iv
	eststo C_2_z_`var': quietly reghdfe z_`var' z_asset_pca  `indvars' if abs(bw)<=10 & group!=., absorb( $fe ) cluster( $clust )

	// D: Changing control districts
	cap drop treat_1
	gen treat_1 = treat_urban
	label variable treat_1 "\hspace{0.1cm} \textit{Reform}"
	eststo D_1_`var': capture quietly reghdfe `var' b_cert $ctrl if abs(bw)<=10 & g_urban!=., absorb( $fe ) cluster( $clust )	
	quietly ivreghdfe `var' (b_cert=treat_1) $ctrl i.yob if abs(bw)<=10 & g_urban!=., a( $clust ) cluster( $clust )
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		matrix define bpmat = J(1,1,0)
		matrix colnames bpmat = "b_cert"
		boottest b_cert, reps(5000) seed(94305) noci statistic(t)
		matrix bpmat[1,1] = `r(p)'
		estadd matrix pboot = bpmat
		est sto D_1_`var'_iv
	eststo D_1_z_`var': quietly reghdfe z_`var' z_asset_pca $ctrl if abs(bw)<=10 & g_urban!=., absorb( $fe ) cluster( $clust )
		
	cap drop treat_1
	gen treat_1 = treat_all
	label variable treat_1 "\hspace{0.1cm} \textit{Reform}"	
	eststo D_2_`var': capture quietly reghdfe `var' b_cert $ctrl if abs(bw)<=10, absorb( $fe ) cluster( $clust )	
	quietly ivreghdfe `var' (b_cert=treat_1) $ctrl i.yob if abs(bw)<=10, a( $clust ) cluster( $clust )
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		matrix define bpmat = J(1,1,0)
		matrix colnames bpmat = "b_cert"
		boottest b_cert, reps(5000) seed(94305) noci statistic(t)
		matrix bpmat[1,1] = `r(p)'
		estadd matrix pboot = bpmat
		est sto D_2_`var'_iv
	eststo D_2_z_`var': quietly reghdfe z_`var' z_asset_pca $ctrl if abs(bw)<=10, absorb( $fe ) cluster( $clust )
		
	// E: Varying fixed effects	
	eststo E_1_`var': capture quietly reghdfe `var' b_cert $ctrl if abs(bw)<=10 & group!=., absorb( yob dcode_2012 ) cluster( dcode_2012 )	
	quietly ivreghdfe `var' (b_cert=treat) $ctrl i.yob if abs(bw)<=10 & group!=., a( dcode_2012 ) cluster( dcode_2012 )
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		matrix define bpmat = J(1,1,0)
		matrix colnames bpmat = "b_cert"
		boottest b_cert, reps(5000) seed(94305) noci statistic(t)
		matrix bpmat[1,1] = `r(p)'
		estadd matrix pboot = bpmat
		est sto E_1_`var'_iv
	eststo E_1_z_`var': quietly reghdfe z_`var' z_asset_pca $ctrl if abs(bw)<=10 & group!=., absorb( yob dcode_2012 ) cluster( dcode_2012 )

	eststo E_2_`var': capture quietly reghdfe `var' b_cert $ctrl if abs(bw)<=10 & group!=., absorb( $fe ) cluster( dcode_1967_y )	
	quietly ivreghdfe `var' (b_cert=treat) $ctrl i.yob if abs(bw)<=10 & group!=., a( $clust ) cluster( dcode_1967_y )
		sum `var' if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
		matrix define bpmat = J(1,1,0)
		matrix colnames bpmat = "b_cert"
		boottest b_cert, reps(5000) seed(94305) noci statistic(t)
		matrix bpmat[1,1] = `r(p)'
		estadd matrix pboot = bpmat
		est sto E_2_`var'_iv		
	eststo E_2_z_`var': quietly reghdfe z_`var' z_asset_pca $ctrl if abs(bw)<=10 & group!=., absorb( $fe ) cluster( dcode_1967_y )
		
}

foreach spec in A B C D E {
	
	foreach v in 1 2 {
		
		local titles "& \multicolumn{3}{c}{I. Education} & \multicolumn{3}{c}{II. Social security} \\"
		local lines "\cmidrule(lr{0.5em}){2-4} \cmidrule(lr{0.5em}){5-7} "
		local names " & Pri. & Sec. & Uni. & HI & Priv. & State \\"
		local numbers " & (1) & (2) & (3) & (4) & (5) & (6) \\ \midrule"

		esttab `spec'_`v'_ed_primary `spec'_`v'_ed_secondary `spec'_`v'_ed_uni `spec'_`v'_s_nih `spec'_`v'_s_pens_priv `spec'_`v'_s_pens_pub using "Additional results/Table 2_`spec'`v'.tex", /// 
			coeflabel(b_cert "\hspace{0.1cm} Registered ($\beta^{OLS})$") replace nonumber posthead("`titles'" "`lines'" "`names'" "`numbers'") nomtitles b(2) se label f nogaps booktabs sfmt(0) /// 
			noobs compress nobaselevels keep(b_cert) star(* 0.10 ** 0.05 *** 0.01)	
		esttab `spec'_`v'_ed_primary_iv `spec'_`v'_ed_secondary_iv `spec'_`v'_ed_uni_iv `spec'_`v'_s_nih_iv `spec'_`v'_s_pens_priv_iv `spec'_`v'_s_pens_pub_iv using "Additional results/Table 2_`spec'`v'.tex", append cells(b(star pvalue(pboot) fmt(%9.2f)) se(par fmt(%9.2f)) pboot(par([ ]) fmt(%9.2f))) coeflabel(b_cert  "\hspace{0.1cm} $\widehat{\text{Registered}}$ ($\beta^{IV})$") nonumber nomtitles label f nogaps booktabs noobs stats(DV_Mean widstat, labels("\midrule DV Mean" "FS F-statistic") layout("@") fmt(2 1)) compress nobaselevels keep(b_cert) nolines prehead(\midrule) star(* 0.10 ** 0.05 *** 0.01) mlabels(,none) collabels(,none)
		esttab `spec'_`v'_z_ed_primary `spec'_`v'_z_ed_secondary `spec'_`v'_z_ed_uni `spec'_`v'_z_s_nih `spec'_`v'_z_s_pens_priv `spec'_`v'_z_s_pens_pub using "Additional results/Table 2_`spec'`v'.tex", /// 
			coeflabel(z_asset_pca "\$\rho\$(Wealth, DV)") append nonumber nomtitles b(2) not label f nogaps booktabs sfmt(0) noobs scalars("N Observations") compress nostar nobaselevels keep(z_asset_pca)
			
	}
}

*** Table 4 using every estimating permutation from Table A4 
********************************************************************************

foreach var of varlist jm_* lee_* sch_* {
	cap drop z_urban_`var'
	cap drop z_all_`var'
	egen z_urban_`var' = std(`var') if abs(bw) <= 10 & g_urban!=.
	egen z_all_`var' = std(`var') if abs(bw) <= 10
}

global ctrl = "c.pl_male c.pl_male#(i.dcode_1967_g)"
global het  = "c.z_jm_gdp_cap c.z_jm_pop_dens"

local spec1 "c.treat"
local spec2 "c.treat##(c.z_jm_rate_payers_cap) c.treat#($het) "
local spec3 "c.treat##(c.z_jm_tax_cap) c.treat#($het) "
local spec4 "c.treat##(c.z_lee_diff_tax) c.treat#($het) "
local spec5 "c.treat##(c.z_jm_rate_payers_cap c.z_jm_tax_cap c.z_lee_diff_tax) c.treat#($het) "
local spec6 "c.treat##(c.z_sch_primary_n_66) c.treat#($het) "
local spec7 "c.treat##(c.z_sch_secondary_any_66) c.treat#($het) "
local spec8 "c.treat##(c.z_sch_primary_n_66 c.z_sch_secondary_any_66) c.treat#($het) "

label variable z_urban_jm_rate_payers_cap "Share paying tax"
label variable z_urban_jm_tax_cap "Tax per capita"
label variable z_urban_lee_diff_tax "\$(\tau^{\text{Max}} - \tau^{\text{Min}})\$"
label variable z_urban_sch_primary_n_66 "Primary schools"
label variable z_urban_sch_secondary_any_66 "Secondary school"

label variable z_all_jm_rate_payers_cap "Share paying tax"
label variable z_all_jm_tax_cap "Tax per capita"
label variable z_all_lee_diff_tax "\$(\tau^{\text{Max}} - \tau^{\text{Min}})\$"
label variable z_all_sch_primary_n_66 "Primary schools"
label variable z_all_sch_secondary_any_66 "Secondary school"

label variable z_jm_rate_payers_cap "Share paying tax"
label variable z_jm_tax_cap "Tax per capita"
label variable z_lee_diff_tax "\$(\tau^{\text{Max}} - \tau^{\text{Min}})\$"
label variable z_sch_primary_n_66 "Primary schools"
label variable z_sch_secondary_any_66 "Secondary school"
label variable z_jm_gdp_cap "GDP per capita"
label variable z_jm_pop_dens "Population density"

global het_urban  = "c.z_urban_jm_gdp_cap c.z_urban_jm_pop_dens"
local spec1D1 "c.treat_urban"
local spec2D1 "c.treat_urban##(c.z_urban_jm_rate_payers_cap) c.treat_urban#($het_urban) "
local spec3D1 "c.treat_urban##(c.z_urban_jm_tax_cap) c.treat_urban#($het_urban) "
local spec4D1 "c.treat_urban##(c.z_urban_lee_diff_tax) c.treat_urban#($het_urban) "
local spec5D1 "c.treat_urban##(c.z_urban_jm_rate_payers_cap c.z_urban_jm_tax_cap c.z_urban_lee_diff_tax) c.treat_urban#($het_urban) "
local spec6D1 "c.treat_urban##(c.z_urban_sch_primary_n_66) c.treat_urban#($het_urban) "
local spec7D1 "c.treat_urban##(c.z_urban_sch_secondary_any_66) c.treat_urban#($het_urban) "
local spec8D1 "c.treat_urban##(c.z_urban_sch_primary_n_66 c.z_urban_sch_secondary_any_66) c.treat_urban#($het_urban) "

global het_all  = "c.z_all_jm_gdp_cap c.z_all_jm_pop_dens"
local spec1D2 "c.treat_all"
local spec2D2 "c.treat_all##(c.z_all_jm_rate_payers_cap) c.treat_all#($het_all) "
local spec3D2 "c.treat_all##(c.z_all_jm_tax_cap) c.treat_all#($het_all) "
local spec4D2 "c.treat_all##(c.z_all_lee_diff_tax) c.treat_all#($het_all) "
local spec5D2 "c.treat_all##(c.z_all_jm_rate_payers_cap c.z_all_jm_tax_cap c.z_all_lee_diff_tax) c.treat_all#($het_all) "
local spec6D2 "c.treat_all##(c.z_all_sch_primary_n_66) c.treat_all#($het_all) "
local spec7D2 "c.treat_all##(c.z_all_sch_secondary_any_66) c.treat_all#($het_all) "
local spec8D2 "c.treat_all##(c.z_all_sch_primary_n_66 c.z_all_sch_secondary_any_66) c.treat_all#($het_all) "

foreach v in 1 2 3 4 5 6 7 8 {

	eststo A_1_`v': quietly reghdfe b_cert $ctrl `spec`v'' if abs(bw) <= 5 & group!=., a( $fe ) cluster( $clust )
		sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
	eststo A_2_`v': quietly reghdfe b_cert $ctrl `spec`v'' if group!=., a( $fe ) cluster( $clust )
		sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'		

	eststo B_1_`v': quietly reghdfe b_cert $ctrl `spec`v'' if abs(bw)<=10 & bw!=0 & group!=., a( $fe ) cluster( $clust )
		sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'		
	eststo B_2_`v': quietly reghdfe b_cert $ctrl `spec`v'' if abs(bw)<=10 & age!=40&age!=45&age!=50&age!=55 & group!=., a( $fe ) cluster( $clust )
		sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'		
	
	local distvars "c.z_jm_s_enrollment_cap c.z_jm_hh_size c.z_jm_h_beds_cap c.z_jm_pop_growth c.z_jm_emp_ag_share"	
	eststo C_1_`v': quietly reghdfe b_cert $ctrl c.post#(`distvars') `spec`v'' if abs(bw)<=10 & group!=., a( $fe ) cluster( $clust )
		sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'			
	local indvars "c.pl_* i.dcode_2012#(c.pl_*)"
	eststo C_2_`v': quietly reghdfe b_cert `indvars'  `spec`v'' if abs(bw)<=10 & group!=., a( $fe ) cluster( $clust )
		sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'			
		
	eststo D_1_`v': quietly reghdfe b_cert $ctrl `spec`v'D1' if abs(bw)<=10 & g_urban!=., a( $fe ) cluster( $clust )
		sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'		
	eststo D_2_`v': quietly reghdfe b_cert $ctrl `spec`v'D2' if abs(bw)<=10, a( $fe ) cluster( $clust )
		sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'	
		
	eststo E_1_`v': quietly reghdfe b_cert $ctrl `spec`v'' if abs(bw)<=10 & group!=., a( dcode_2012 yob ) cluster( dcode_2012 )
		sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'			
	eststo E_2_`v': quietly reghdfe b_cert $ctrl `spec`v'' if abs(bw)<=10 & group!=., a( $fe ) cluster( dcode_1967_y )
		sum b_cert if e(sample)==1
		estadd scalar DV_Mean = `r(mean)'
}

cd "$output"

foreach spec in A B C D E {
	
	foreach v in 1 2 {
		
		local over "& & \multicolumn{4}{c}{I. Taxation} & \multicolumn{3}{c}{II. Local public goods} \\"
		local lines "\cmidrule(lr{0.5em}){3-6} \cmidrule(lr{0.5em}){7-9}"
		local numbers "& (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8) \\"

		esttab `spec'_`v'_1 `spec'_`v'_2 `spec'_`v'_3 `spec'_`v'_4 `spec'_`v'_5 `spec'_`v'_6 `spec'_`v'_7 `spec'_`v'_8 using "Additional results/Table 4_`spec'`v'.tex", ///
			replace posthead("`over'" "`lines'" "`numbers'") nomtitles b(3) se label f nogaps booktabs scalars("DV_Mean Outcome mean" "N Observations") /// 
			sfmt(2 0) noobs nonumber compress star(* 0.10 ** 0.05 *** 0.01) nobaselevels  substitute(\_ _) /// 
			keep(trea* *trea*#c.z_*) drop(*trea*#c.z_*gdp_* *trea*#c.z_*pop_*)
	}
}
